#######################################
# Replication file for Chaco Province #
#######################################


rm(list=ls(all=TRUE))
#setwd("")
library(eiPack)
library(coda)
library(foreign)

data=read.csv(file.path(getwd(),"data","CHACO_data.csv"))

# Define row and column marginals of the transition matrix.
GFPV=data$gral_501
GUPC=data$gral_519
GPO=data$gral_14
GN=data$votantes-GFPV-GUPC-GPO


PFPVW=data$paso_3118
PFPVL=data$paso_3122+data$paso_3120+data$paso_3119
PUPCA=data$paso_3123
PUPCB=data$paso_3124
PPO=data$paso_14
POther=data$paso_60+data$paso_182+data$paso_518+data$paso_179
PNV=data$votantes-PFPVW-PFPVL-PUPCA-PUPCB-PPO-POther

data1=data.frame(PFPVW,PFPVL,PUPCA,PUPCB,PPO,POther,PNV,GFPV,GUPC,GPO,GN)



# Tunning the EI algorithm
tune.nocov=tuneMD(cbind(GFPV,GUPC,GPO,GN)~cbind(PFPVW,PFPVL,PUPCA,PUPCB,PPO,POther,PNV),
                  data=data1,ntunes=1,totaldraws=200000)



# MCMC
out.nocov=ei.MD.bayes(cbind(GFPV,GUPC,GPO,GN)~cbind(PFPVW,PFPVL,PUPCA,PUPCB,PPO,POther,PNV),
                      covariate=NULL,data=data1,tune.list=tune.nocov,verbose=0,ret.beta="r",thin=80,burnin=100000,ret.mcmc=TRUE) 



# Obtain transition matrix for each mesa: 
source(file.path(getwd(),"Codes","getbetas.R"))
beta11=getbetas(out.nocov,"PFPVW","GFPV")
beta12=getbetas(out.nocov,"PFPVW","GUPC")
beta13=getbetas(out.nocov,"PFPVW","GPO")
beta14=getbetas(out.nocov,"PFPVW","GN")
beta21=getbetas(out.nocov,"PFPVL","GFPV")
beta22=getbetas(out.nocov,"PFPVL","GUPC")
beta23=getbetas(out.nocov,"PFPVL","GPO")
beta24=getbetas(out.nocov,"PFPVL","GN")
beta31=getbetas(out.nocov,"PUPCA","GFPV")
beta32=getbetas(out.nocov,"PUPCA","GUPC")
beta33=getbetas(out.nocov,"PUPCA","GPO")
beta34=getbetas(out.nocov,"PUPCA","GN")
beta41=getbetas(out.nocov,"PUPCB","GFPV")
beta42=getbetas(out.nocov,"PUPCB","GUPC")
beta43=getbetas(out.nocov,"PUPCB","GPO")
beta44=getbetas(out.nocov,"PUPCB","GN")
beta51=getbetas(out.nocov,"PPO","GFPV")
beta52=getbetas(out.nocov,"PPO","GUPC")
beta53=getbetas(out.nocov,"PPO","GPO")
beta54=getbetas(out.nocov,"PPO","GN")
beta61=getbetas(out.nocov,"POther","GFPV")
beta62=getbetas(out.nocov,"POther","GUPC")
beta63=getbetas(out.nocov,"POther","GPO")
beta64=getbetas(out.nocov,"POther","GN")
beta71=getbetas(out.nocov,"PNV","GFPV")
beta72=getbetas(out.nocov,"PNV","GUPC")
beta73=getbetas(out.nocov,"PNV","GPO")
beta74=getbetas(out.nocov,"PNV","GN")

betaout=data.frame(beta11,beta12,beta13,beta14,
                   beta21,beta22,beta23,beta24,
                   beta31,beta32,beta33,beta34,
                   beta41,beta42,beta43,beta44,
                   beta51,beta52,beta53,beta54,
                   beta61,beta62,beta63,beta64,
                   beta71,beta72,beta73,beta74,
                   data$votantes,data$prov,data$comuna,data$circ,data$mesa)

write.table(betaout,file=file.path(getwd(),"Output","CHACO_Mesa.txt"),col.names=NA)

# Obtain aggregate transition matrix. 
source(file.path(getwd(),"Codes","getbetasagr.R"))
beta11=getbetasagr(out.nocov,"PFPVW","GFPV")
beta12=getbetasagr(out.nocov,"PFPVW","GUPC")
beta13=getbetasagr(out.nocov,"PFPVW","GPO")
beta14=getbetasagr(out.nocov,"PFPVW","GN")
beta21=getbetasagr(out.nocov,"PFPVL","GFPV")
beta22=getbetasagr(out.nocov,"PFPVL","GUPC")
beta23=getbetasagr(out.nocov,"PFPVL","GPO")
beta24=getbetasagr(out.nocov,"PFPVL","GN")
beta31=getbetasagr(out.nocov,"PUPCA","GFPV")
beta32=getbetasagr(out.nocov,"PUPCA","GUPC")
beta33=getbetasagr(out.nocov,"PUPCA","GPO")
beta34=getbetasagr(out.nocov,"PUPCA","GN")
beta41=getbetasagr(out.nocov,"PUPCB","GFPV")
beta42=getbetasagr(out.nocov,"PUPCB","GUPC")
beta43=getbetasagr(out.nocov,"PUPCB","GPO")
beta44=getbetasagr(out.nocov,"PUPCB","GN")
beta51=getbetasagr(out.nocov,"PPO","GFPV")
beta52=getbetasagr(out.nocov,"PPO","GUPC")
beta53=getbetasagr(out.nocov,"PPO","GPO")
beta54=getbetasagr(out.nocov,"PPO","GN")
beta61=getbetasagr(out.nocov,"POther","GFPV")
beta62=getbetasagr(out.nocov,"POther","GUPC")
beta63=getbetasagr(out.nocov,"POther","GPO")
beta64=getbetasagr(out.nocov,"POther","GN")
beta71=getbetasagr(out.nocov,"PNV","GFPV")
beta72=getbetasagr(out.nocov,"PNV","GUPC")
beta73=getbetasagr(out.nocov,"PNV","GPO")
beta74=getbetasagr(out.nocov,"PNV","GN")

betaoutagr=data.frame(beta11,beta12,beta13,beta14,
                   beta21,beta22,beta23,beta24,
                   beta31,beta32,beta33,beta34,
                   beta41,beta42,beta43,beta44,
                   beta51,beta52,beta53,beta54,
                   beta61,beta62,beta63,beta64,
                   beta71,beta72,beta73,beta74,
                   data$votantes,data$prov,data$comuna,data$circ,data$mesa)

betaoutagr=as.vector(betaoutagr[1,1:100])
betaoutagr=matrix(betaoutagr,4)
betaoutagr_coefs=t(matrix(betaoutagr[1,],5))
betaoutagr_sd=t(matrix(betaoutagr[2,],5))
write.table(betaoutagr_coefs,file=file.path(getwd(),"Output","CHACOAggregate_Coefs.txt"),col.names=NA)
write.table(betaoutagr_sd,file=file.path(getwd(),"Output","CHACOAggregate_SD.txt"),col.names=NA)


#save(out.nocov, file = "CHACOMCMC.Rdata") # Saves MCMC runs


#Geweke Convergence Statistics 
geweke=geweke.diag(out.nocov$draws$Beta)
gewekestats=geweke[[1]]
gewekedf=data.frame(gewekestats)
#write.table(gewekedf,file="CHACOgewekelist.txt",col.names=NA)


jpeg(file.path(getwd(),"Output","CHACOHistogramGeweke.jpeg"))
hist(gewekestats)
dev.off()

print("Max Geweke")
print(max(gewekestats))
print("Min Geweke")
print(min(gewekestats))
print("Percentage Above 2.35")
print(length(gewekestats[abs(gewekestats)>2.35])/length(gewekestats))
print("Percentage Above 1.96")
print(length(gewekestats[abs(gewekestats)>1.96])/length(gewekestats))
print("Percentage Above 1.5")
print(length(gewekestats[abs(gewekestats)>1.5])/length(gewekestats))
print("Percentage Above 1.5")
print(length(gewekestats[abs(gewekestats)>1])/length(gewekestats))
print("Percentage Above 0.5")
print(length(gewekestats[abs(gewekestats)>0.5])/length(gewekestats))
